#include "matrixMul.h"
#include <cuda_runtime.h>
#include <stdlib.h>
#include <stdio.h>

void matrixMulCPU(const Matrix A, const Matrix B, Matrix C)
{
    for (unsigned int i = 0; i < C.height; ++i)
        for (unsigned int j = 0; j < C.width; ++j) {
            double sum = 0;
            for (unsigned int k = 0; k < A.width; ++k) {
                double a = A.element[i * A.width + k];
                double b = B.element[k * B.width + j];
                sum += a * b;
            }
            C.element[i * C.width + j] = (float)sum;
        }
}

void randomInitMatrix(Matrix &M)
{
    M.element=(float*)malloc(M.height*M.width*sizeof(float));
    for (unsigned int i=0; i< M.height; i++){
        for (unsigned int j=0; j<M.width; j++){
            M.element[i*M.width+ j ] = rand()/(float) RAND_MAX;
        }
    }
}

void OnesMatrix(Matrix &M)
{
    //M.element=(float*)malloc(M.height*M.width*sizeof(float));
    for (unsigned int i=0; i< M.height; i++){
        for (unsigned int j=0; j<M.width; j++){
            M.element[i*M.width+ j ] = 1;
        }
    }

}

__device__ float GetElement(const Matrix A, int row, int col)
{
    return A.element[row*A.stride+col];
}

__device__ void SetElement(const Matrix A, int row, int col, float value)
{
    A.element[row*A.stride+col]=value;
}

__device__ Matrix GetSubMatrix(Matrix A, int row, int col)
{
    Matrix Asub;
    Asub.width=BLOCK_SIZE;
    Asub.height=BLOCK_SIZE;
    Asub.stride=A.stride;
    Asub.element=&A.element[A.stride*BLOCK_SIZE*row+BLOCK_SIZE*col];

    return Asub;
}

__global__ void matrixMulKernel(Matrix A, Matrix B, Matrix C)
{

    float Cvalue = 0;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    for (int e = 0; e < A.width; ++e){

        Cvalue += A.element[row * A.width + e]
        * B.element[e * B.width + col];
        C.element[row * C.width + col] = Cvalue;

    }


}

void matrixMul(const Matrix A, const Matrix B, Matrix C)
{
    Matrix d_A;
    Matrix d_B;
    Matrix d_C;
    d_A.width=d_A.stride=A.width;
    d_A.height=A.height;
    d_B.width=d_B.stride=B.width;
    d_B.height=B.height;
    d_C.width=d_C.stride=C.width;
    d_C.height=C.height;

    size_t sizeA = A.width*A.height*sizeof(float);
    cudaMalloc((void**)&d_A.element, sizeA);
    cudaMemcpy(d_A.element, A.element, sizeA, cudaMemcpyHostToDevice);

    size_t sizeB = B.width*B.height*sizeof(float);
    cudaMalloc((void**)&d_B.element, sizeB);
    cudaMemcpy(d_B.element, B.element, sizeB, cudaMemcpyHostToDevice);

    size_t sizeC = C.width*C.height*sizeof(float);
    cudaMalloc((void**)&d_C.element, sizeC);
    cudaMemcpy(d_C.element, C.element, sizeC, cudaMemcpyHostToDevice);

    dim3 dimBLOCK(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid(B.width/dimBLOCK.x, A.height/dimBLOCK.y);
    matrixMulKernel<<<dimGrid, dimBLOCK>>>(d_A, d_B, d_C);

    cudaMemcpy(C.element, d_C.element, sizeC, cudaMemcpyDeviceToHost);

    cudaFree(d_A.element);
    cudaFree(d_B.element);
    cudaFree(d_C.element);

}

bool initCUDA()
{
    int count;

    cudaGetDeviceCount(&count);
    if(count == 0) {
        fprintf(stderr, "There is no device.\n");
        return false;
    }

    int i;
    for(i = 0; i < count; i++) {
        cudaDeviceProp prop;
        if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
            if(prop.major >= 1) {
              break;
            }
        }
    }

    if(i == count) {
        fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
        return false;
    }

    cudaSetDevice(i);

    return true;
}

int main()
{
    //initCUDA();

    Matrix A;
    Matrix B;
    Matrix Cgpu;
    Matrix Ccpu;
    A.height = HA;
    A.width  = WA;
    A.stride = WA;
    B.height = HB;
    B.width  = WB;
    B.stride  = WB;
    Ccpu.height = HC;
    Ccpu.width  = WC;
    Ccpu.stride = WC;
    Cgpu.height = HC;
    Cgpu.width  = WC;
    Cgpu.stride = WC;

    randomInitMatrix(A);
    randomInitMatrix(B);
    randomInitMatrix(Ccpu);
    randomInitMatrix(Cgpu);

    //OnesMatrix(A);
    //OnesMatrix(B);

    /*
    randomInitMatrix(A.element, A.height, A.width);
    randomInitMatrix(B.element, B.height, B.width);
    randomInitMatrix(Ccpu.element, Ccpu.height, Ccpu.width);
    randomInitMatrix(Cgpu.element, Cgpu.height, Cgpu.width);
    */

    matrixMulCPU(A, B, Ccpu);
    matrixMul(A,B,Cgpu);

    float max = 0;
    for (unsigned int i=0; i < HC; i++){
        for (unsigned int j=0; j< WC; j++){
            //printf("row:%d col$:%d cpuvalue:%f gpuvalue:%f\n",i,j, Ccpu.element[i*WC+j], Cgpu.element[i*WC+j]);

            float diff = abs(Ccpu.element[i*WC+j]-Cgpu.element[i*WC+j]);
            if (diff > max){
                max = diff;
            }
        }
    }

    free(A.element);
    free(B.element);
    free(Ccpu.element);
    free(Cgpu.element);


    printf("%f\n", max);
}
